function  kappap_matrix = resolving_kappap( p_2D_new,  par  )
% Usage: adjust the eta for all the solution files of Smolayak grid.

%% *** Adjust the kappap function again. 
N_w = par.N_w;    N_lambda = par.N_lambda;
AH=par.AH;   AL=par.AL;   rho=par.rho;  sigmaK = par.sigmaK;   investment=par.investment;   
alpha=par.alpha;   beta=par.beta;
max_jump=par.max_jump;   
w_grid = par.w_grid;   lambda_grid = par.lambda_grid;

p_2D_fitted =  scatteredInterpolant(   kron(ones(N_lambda,1), w_grid ),    kron(lambda_grid, ones(N_w,1)) ,  reshape( p_2D_new, N_w*N_lambda, 1 ),   'linear');
par.p_2D_fitted = p_2D_fitted;

% Other output variables
init_matrix = zeros(  N_w,  N_lambda );
kappap_matrix = init_matrix;  

% error control
epsilon=0.001; % maximum error in solving kappap

for( counter_lambda = 1:N_lambda)
    lambda = lambda_grid(counter_lambda);   par.lambda=lambda;  
%     disp( [   'Inside 2D iteration: lambda=', num2str(lambda)   ]   )
    
    kappap_sol = zeros( N_w, 1 ); xK_sol = zeros(N_w,1);  yK_sol = zeros(N_w,1);    psi_new = zeros(N_w,1);   solved_vec = zeros(N_w,1);
    p_new = p_2D_new(:,counter_lambda) ;      pw_new=derivative( w_grid, p_new );
    %== Step 1: For each grid point of w, solve kappa, xK, yK. 
    for(k = 2:(N_w-1))
        % Decide whether at p, and w we have psi=1
        w=w_grid(k);    p=p_new(k);   pw=pw_new(k);
        yK_fun = @(xK) ( 1 - w*xK ) ./ (1-w);

        sigmap_fun = @(xK) pw.*w.*(1-w).*(xK- yK_fun(xK) ) ./ ( 1- pw.*w.*(1-w).*(xK- yK_fun(xK) )  ).* sigmaK; 
        sigmah_fun = @(xK) yK_fun(xK) .* (sigmaK + sigmap_fun(xK) ); 
        sigma_fun = @(xK)  xK .* ( sigmaK + sigmap_fun(xK) );

        par.yK_fun=yK_fun;  par.sigmap_fun=sigmap_fun;  par.sigmah_fun=sigmah_fun;  par.sigma_fun=sigma_fun;
        par.w_ind = w;    par.p_ind = p;   par.pw_ind = pw;   

        xK_min = 1  ;     % The value of xK such that yK=xK
        if(pw<=0)
            xK_max1 =  1000 ;   % Don't limit the value of xK, and use a uniform set of first order conditions.
        else
            xK_max1 = 1/( pw*w) + 1 ;    % maximum value restricted by positive volatility
        end
        xK_max1 = min( xK_max1,  (1+alpha*beta)/(alpha*beta) ) ;  % Maximum likelihood restricted by kappa equation
        xKs = exp(linspace( log(xK_min+1), log(xK_max1*0.98+1) ,100)') - 1;
        N_sample = 15;
        % kappaps = [0;  exp(linspace( log(0.001), log(0.9), N_sample-1 )')] ; 
        kappaps = linspace( 0, min( max_jump, p-0.02 ) , N_sample )' ; 
        xK_roots = zeros( numel(kappaps), 1 );
        solved_inner = (ones( N_sample,1 )==1);   
        for(i = 1:N_sample)
            kappap = kappaps(i);
            xK_max = min( xK_max1,   (1+alpha*beta)/(kappap+alpha*beta)  );    par.xK_max = xK_max;   % Note: we should use min(xK_max, 1/kappap) so that it updates for each kappap, instead of min(par.xK_max, 1/kappap)
            par.xK_min = xK_min;  
            [ xK_sol , fval ,  ~,  ~] = fsolve( @(xK) 10*xK_iterative_residual(xK, kappap, par) , xK_min*1.1 ,  optimset('Display','off') );  % note: scaling by 10 improves solution accuracy                
            if( abs(fval)>epsilon  )
                xK_corner = 1/w;
                if( xK_iterative_residual(1, kappap, par) * xK_iterative_residual(1.001, kappap, par) < 0 )
                    xK_sol = 1; solved_inner(i)=true;
                elseif( xK_iterative_residual(xK_corner, kappap, par) > 0  &&  xK_corner<=xK_max )   % Corner case works
                    xK_sol = xK_corner; solved_inner(i)=true;
                else   % Corner case does not work
                    disp( [ 'In solving xK1, a solution fails (should not happen)!!! At kappap=', num2str(kappap), '(iteration=', num2str(i)  , '), w=', num2str(w),  '(main iteration=', num2str(k), ')' ] );
                    xK_sol = xK_min;   solved_inner(i)=false;
                end
            end
            xK_roots(i)  =  xK_sol;
        end
        close all;
        % roots_SLM = slmengine( kappaps(solved) ,xK_roots(solved), 'plot', 'on', 'decreasing', 'on',    'knots',   floor( N_sample /2)   );  
        roots_xK = interp1( kappaps(solved_inner) ,xK_roots(solved_inner),'linear','pp');

        % Then add parameters to par for solving kappap
        par.xK_roots_fitted = roots_xK;

        counter=0;   
        kp0=0;       max_kp =  max(kappaps)-0.01   ;  step_size=0.1;    maxit = floor(max_kp/step_size) ;  fval=1;
        while( abs(fval)>epsilon & counter<maxit & kp0<max_kp )   % exitflat<0 means no solution has been found. If so, we increase the initial value kp0 and search again.
                counter = counter+1;
                [kappap_sol(k), fval ,~,~] = fsolve( @(kp) kappa_iterative_residual(kp, par),  kp0, optimset( 'Display','none' ) );  % Important: use subscript k
                kp0=kp0+step_size;  % Note: it is very important to progressively increase the initial point. If we plot the residual as a function of kappa, then we can observe there are multiple solutions. We pick the solution closest to zero by this algorithm. 
        end
        solved= (abs(fval)<epsilon) ;
        if( ~solved )  
            disp( [ 'In solving kappap, a solution fails! At w=', num2str(w),  '(iteration=', num2str(k), ')' ] )
        end
        solved_vec(k) = solved;
        xK_sol(k) = ppval( roots_xK,kappap_sol(k));  yK_sol(k) = yK_fun( xK_sol(k)  );
        psi_new(k) = w.*xK_sol(k) ; 
    end
    %== Step 2: Get updated results for price function. Note that psi<=1 must be satisfied.
    psi_sol = min( psi_new, 1 ) ;        psi_sol( w_grid==1 ) = 1;    psi_sol( w_grid==0 ) = 0; 
    p_sol_next = zeros(  numel(w_grid) ,1 );  
    for(i = 1:numel(w_grid))
        p_sol_next(i) = fsolve(  @(p)  investment(p) + rho*p - AL - psi_sol(i)*( AH-AL ),  1 ,   optimoptions( 'fsolve','Display','none' ) );  
    end
    % We need to smooth the function around zero
    solved_vec=(solved_vec==1 & (~isnan(solved_vec)) );  solved_vec(1)=true;  solved_vec(N_w)=true;
 
    % Update p_vecplot( lambda_grid, kappap_matrix )

    kappap_vec = interp1(w_grid(solved_vec) ,  kappap_sol(solved_vec),  w_grid);
    kappap_matrix(:,counter_lambda) = kappap_vec;  
end

